Glasser regions list
glasser.regions <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist.csv")
glasser.snr.exclude <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/SNR/glasser_SNR_exclusion.csv")
glasser.frontal <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist_frontallobe.csv")
glasser.lh.frontallobe.order <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/LeftDistmat_frontallobe_regionorder.txt", header = F) %>% set_names("orig_parcelname")
glasser.rh.frontallobe.order <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/RightDistmat_frontallobe_regionorder.txt", header = F) %>% set_names("orig_parcelname")S-A Axis
SAaxis.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii")
SAaxis <- data.frame(SA.axis = rank(SAaxis.cifti$data), orig_parcelname = names(SAaxis.cifti$Parcel))BigBrain cytoarchitectural variation
bigbrain.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/BigBrain_histologygradient/BigBrain.Histology.pscalar.nii")
bigbrain <- data.frame(cytoarchitecture.gradient = bigbrain.cifti$data, orig_parcelname = names(bigbrain.cifti$Parcel))Myelin development data for correlation significance testing
regional.rate <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_depth_development/regional_averagederivative_frontallobe.csv")
regional.maturationage <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_depth_development/regional_maturationage_frontallobe.csv")Autocorrelation-preserving surrogates
#Hemisphere-specific surrogates generated by brainsmash (n = 500 per hemisphere)
surrogate_maps_rh <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/SAaxis.rh.frontallobe.surrogates.txt") %>% t() %>% as.data.frame() #72 rows are RH regions, 500 columns are 500 surrogate maps
surrogate_maps_lh <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/SAaxis.lh.frontallobe.surrogates.txt") %>% t() %>% as.data.frame() #72 rows are LH regions, 500 columns are 500 surrogate maps
#Reflect data across hemisphere to generate frontal lobe surrogate maps
surrogate_maps_rh <- rbind(surrogate_maps_rh, surrogate_maps_rh)
surrogate_maps_lh <- rbind(surrogate_maps_lh, surrogate_maps_lh)
#Add in region names in correct order
glasser.frontallobe.order <- rbind(glasser.rh.frontallobe.order, glasser.lh.frontallobe.order)
surrogate_maps_rh$orig_parcelname <- glasser.frontallobe.order$orig_parcelname
surrogate_maps_lh$orig_parcelname <- glasser.frontallobe.order$orig_parcelname
#NA out low SNR parcels to exclude from null correlation tests
surrogate_maps_rh <- surrogate_maps_rh[!(surrogate_maps_rh$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
surrogate_maps_rh <- left_join(glasser.frontallobe.order, surrogate_maps_rh, by = "orig_parcelname")
surrogate_maps_lh <- surrogate_maps_lh[!(surrogate_maps_lh$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
surrogate_maps_lh <- left_join(glasser.frontallobe.order, surrogate_maps_lh, by = "orig_parcelname")
surrogate_maps <- merge(surrogate_maps_rh, surrogate_maps_lh, by = c("orig_parcelname"), sort = F) #1000 nulls!
colnames(surrogate_maps)[2:1001] <- sprintf("map%s",seq(from = 1, to = 1000))
surrogate_maps <- left_join(surrogate_maps, glasser.regions, by = c("orig_parcelname"))for(map in c("map1", "map100", "map200", "map300", "map400", "map500", "map600", "map700", "map800", "map900", "map100")){
null.map <- surrogate_maps %>% filter(label != "lh_L_10pp") %>% select(label, map)
colnames(null.map)<- c("label", "nullmap")
SA.nullmap.plot <- ggseg(.data = null.map, atlas = "glasser", mapping = aes(fill = nullmap), colour=I("gray50"), size=I(.06)) +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill")
print(SA.nullmap.plot)
}#Function to compute the brainsmash-based permutation p-value (pSMASH!) given an empirical correlation between depth-specific development map and the SA-preserving null S-A axis maps #super smash bros
brainsmash.results <- function(dev.df, dev.metric, this.depth){
df <- dev.df %>% filter(depth == this.depth)
df <- df %>% select(label, any_of(dev.metric))
surrogate_maps.dev <- left_join(surrogate_maps, df, by = c("label")) #ensure surrogate map and dev df regional ordering match
surrogate_maps.dev <- left_join(surrogate_maps.dev, SAaxis, by = c("orig_parcelname"))
surrogate_maps.null <- surrogate_maps.dev %>% select(contains("map")) #1000 null maps
surrogate_maps.empirical <- surrogate_maps.dev %>% select(SA.axis, any_of(dev.metric)) #true data
empirical.cor <- cor(surrogate_maps.empirical$SA.axis, surrogate_maps.empirical[dev.metric], method = c("spearman"), use = "complete.obs") #true correlation
null.cors <- cor(surrogate_maps.empirical[dev.metric], surrogate_maps.null, method = c("spearman"), use = "complete.obs") #correlation with brainsmash based surrogates
mean.nullcor <- mean(null.cors)
if(empirical.cor > 0){
brainsmash.p <- sum(null.cors > empirical.cor[1])/1000
}
if(empirical.cor < 0){
brainsmash.p <- sum(null.cors < empirical.cor[1])/1000
}
smash.results <- list(empirical.cor[1], mean.nullcor, brainsmash.p)
names(smash.results) <- list("empirical.cor", "average.null.cor", "brainsmash.pvalue")
return(smash.results)
}R1 Rate of Increase - S-A Axis Correlation by Depth
## $empirical.cor
## [1] -0.5096381
##
## $average.null.cor
## [1] 0.01730044
##
## $brainsmash.pvalue
## [1] 0.046
## $empirical.cor
## [1] -0.2408838
##
## $average.null.cor
## [1] 0.02411319
##
## $brainsmash.pvalue
## [1] 0.254
Autocorrelation-preserving surrogates
#Hemisphere-specific surrogates generated by brainsmash (n = 500 per hemisphere)
surrogate_maps_rh <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/Bigbrain.rh.frontallobe.surrogates.txt") %>% t() %>% as.data.frame() #72 rows are RH regions, 500 columns are 500 surrogate maps
surrogate_maps_lh <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/Bigbrain.lh.frontallobe.surrogates.txt") %>% t() %>% as.data.frame() #72 rows are LH regions, 500 columns are 500 surrogate maps
#Reflect data across hemisphere to generate frontal lobe surrogate maps
surrogate_maps_rh <- rbind(surrogate_maps_rh, surrogate_maps_rh)
surrogate_maps_lh <- rbind(surrogate_maps_lh, surrogate_maps_lh)
#Add in region names in correct order
glasser.frontallobe.order <- rbind(glasser.rh.frontallobe.order, glasser.lh.frontallobe.order)
surrogate_maps_rh$orig_parcelname <- glasser.frontallobe.order$orig_parcelname
surrogate_maps_lh$orig_parcelname <- glasser.frontallobe.order$orig_parcelname
#NA out low SNR parcels to exclude from null correlation tests
surrogate_maps_rh <- surrogate_maps_rh[!(surrogate_maps_rh$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
surrogate_maps_rh <- left_join(glasser.frontallobe.order, surrogate_maps_rh, by = "orig_parcelname")
surrogate_maps_lh <- surrogate_maps_lh[!(surrogate_maps_lh$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
surrogate_maps_lh <- left_join(glasser.frontallobe.order, surrogate_maps_lh, by = "orig_parcelname")
surrogate_maps <- merge(surrogate_maps_rh, surrogate_maps_lh, by = c("orig_parcelname"), sort = F) #1000 nulls!
colnames(surrogate_maps)[2:1001] <- sprintf("map%s",seq(from = 1, to = 1000))
surrogate_maps <- left_join(surrogate_maps, glasser.regions, by = c("orig_parcelname"))for(map in c("map1", "map100", "map200", "map300", "map400", "map500", "map600", "map700", "map800", "map900", "map100")){
null.map <- surrogate_maps %>% filter(label != "lh_L_10pp") %>% select(label, map)
colnames(null.map)<- c("label", "nullmap")
BB.nullmap.plot <- ggseg(.data = null.map, atlas = "glasser", mapping = aes(fill = nullmap), colour=I("gray50"), size=I(.06)) +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill")
print(BB.nullmap.plot)
}#Function to compute the brainsmash-based permutation p-value (pSMASH!) given an empirical correlation between depth-specific development map and the SA-preserving null BigBrain maps #big brain smash
brainsmash.results <- function(dev.df, dev.metric, this.depth){
df <- dev.df %>% filter(depth == this.depth)
df <- df %>% select(label, any_of(dev.metric))
surrogate_maps.dev <- left_join(surrogate_maps, df, by = c("label")) #ensure surrogate map and dev df regional ordering match
surrogate_maps.dev <- left_join(surrogate_maps.dev, bigbrain, by = c("orig_parcelname"))
surrogate_maps.null <- surrogate_maps.dev %>% select(contains("map")) #1000 null maps
surrogate_maps.empirical <- surrogate_maps.dev %>% select(cytoarchitecture.gradient, any_of(dev.metric)) #true data
empirical.cor <- cor(surrogate_maps.empirical$cytoarchitecture.gradient, surrogate_maps.empirical[dev.metric], method = c("spearman"), use = "complete.obs") #true correlation
null.cors <- cor(surrogate_maps.empirical[dev.metric], surrogate_maps.null, method = c("spearman"), use = "complete.obs") #correlation with brainsmash based surrogates
mean.nullcor <- mean(null.cors)
if(empirical.cor > 0){
brainsmash.p <- sum(null.cors > empirical.cor[1])/1000
}
if(empirical.cor < 0){
brainsmash.p <- sum(null.cors < empirical.cor[1])/1000
}
smash.results <- list(empirical.cor[1], mean.nullcor, brainsmash.p)
names(smash.results) <- list("empirical.cor", "average.null.cor", "brainsmash.pvalue")
return(smash.results)
}R1 Rate of Increase - BigBrain Cytoarchitectural Variation Correlation by Depth
## $empirical.cor
## [1] 0.6371086
##
## $average.null.cor
## [1] -0.004780919
##
## $brainsmash.pvalue
## [1] 0.013
## $empirical.cor
## [1] 0.5435518
##
## $average.null.cor
## [1] -0.004175066
##
## $brainsmash.pvalue
## [1] 0.1